COMPUTING THE REFINED STABILITY CONDITION 
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Abstract. The classical (inviscid) stability analysis of shock waves is based on the Lopatin- 
skh determinant, A — a function of frequencies whose zeros determine the stability of the 
underlying shock. A careful analysis of A shows that in some cases the stable and unstable 
regions of parameter space are separated by an open set of parameters. Zumbrun and Serre 
[Indiana Univ. Math. J., 48 (1999) 937-992] have shown that, by taking account of viscous 
effects not present in the definition of A, it is possible to determine the precise location in 
the open, neutral set of parameter space at which stability is lost. In particular, they show 
that the transition to instability under suitably localized perturbations is determined by an 
"effective viscosity" coefficient. Here, in the simplest possible setting, we propose and im- 
plement two new approaches toward the practical computation of this coefficient. Moreover, 
in a special case, we derive an exact solution of the relevant differential equations. 



1. Introduction 
1.1. Inviscid stability. Consider a planar shock solution, 

u( X ,t) = h> xi>st (i.i) 

\u_, X\ < st 

of a hyperbolic system of conservation laws in d space dimensions: 

d 

d t u + Y,d X] r(u)=0. (1.2) 

3=1 

Here, the unknown u belongs to % C M. n , the state space, and is a function of space 
x = (x\, . . . ,Xd) G K d and time t G i The fluxes / J ' are IR n - valued functions on & . In 
(11. ip . u± are constant states; they are related to the shock speed s by the Rankine-Hugoniot 
condition 

s[u] = [f 1 (u)}, 

where — here and below — square brackets indicate the jump. That is, for any function h of 
the state u, [h(u)] := h(u + ) — h(u_). The linear stability analysis of such solutions is now 
classical; it originates in the studies of pioneers D'yakov [6 J and Erpenbeck [7\. We recall 
that the the centerpiece of the analysis is the Lopatinskii determinant, 

A : {A G C : Re A > 0} x R^ 1 -»■ C , 

a function of frequencies (A, £), where A = 7 + ir e C is dual to time and £ G R d_1 is dual to 
the transverse spatial directions. Zeros of A with 7 = Re A > correspond to perturbations 
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which grow exponentially in time, and, evidently, nonvanishing of A on {7 > 0} is necessary 
for (linearized, inviscid) stability of the shockEI. 

A particularly important example occurs when (11. 2p is the Euler equations of gas dynamics 
and the solution (II. ip represents a planar gas-dynamical shock wave. In this example, the 
vector u would have as its components the mass density, the momentum densities in each of 
the d spatial directions, and the total energy density. In this case, an analysis of A, see [3J or 
the appendix of [13], yields clean, explicit stability criteria in terms of the basic properties 
of the shock and the gaa^ in question. However, Majda [H] has observed that in the setting 
of gas dynamics it sometimes happens that A vanishes on the boundary {7 = 0} but not for 
7 > 0, and in this case we say that the shock is neutrally or weakly stable. Moreover, this 
neutral stability can persist on an open set in parameter space. Physically, the presence of 
neutral zeros of A corresponds to surface or boundary waves, and these zeros are associated 
with a loss of smoothness of the perturbed shock front. Further complicating matters, 
Barmin & Egorushkin [2] have pointed out that experimentally observed instabilities of gas- 
dynamical shocks sometimes occur within the region of weak stability. That is, instabilities 
are sometimes observed in the interior of the weakly stable regime. They postulated that 
the inclusion of nonlinear effects and/or other neglected physical effects would be required 
for the theory to capture the phenomena revealed by these experiments. 

1.2. Viscous stability. By considering the viscous regularization of (II. 2p . 

d d 

d t u + J2d Xj f j (u) Xj =vJ2 d Xj (Bi k (u)d Xk u) , (1.3) 

j=l j,k=l 

Zumbrun & Serre [TH] have shown that the precise location of the transition to instability is 
entirely determined by viscous effects which are neglected in the construction of A. See also 
the more recent paper of Benzoni-Gavage, Serre, & Zumbrun [5]; this latter paper forms the 
foundation of the calculations we present here^. As a point of reference, we note that in the 
aforementioned example of gas dynamics, equation (II. 3p would correspond to the Navier- 
Stokes equations of compressible gas dynamics; the second-order terms on the right-hand 
side of (11.31) model the effects of viscosity and heat conductivity. 

We associate to the solution (II. ip of the inviscid equations (II. 2p a planar viscous profile 
for (11.31) . The viscous profile is a solution of (II. 3p of the form 

u(x,t) = u ( — J , lim u(z) = u ± . (1-4) 

y v ) z->±oo 

Zumbrun & Serre's derivation is based on a low-frequency analysis of the Evans function, 
-D(A, £), associated with ( II .4p . Analogous to the Lopatinskii determinant, the Evans function 
is a spectral determinant; its zero set carries stability information for the planar viscous shock 
wave. In particular, the refined stability condition they derive is given in terms of the sign 

a Indeed, A is homogenous degree one. Thus, any unstable zero generates instabilities of all orders. That 
is, any such instability is of Hadamard type; these instabilities are so violent that these waves will never be 
seen in practice. 

b The stability criteria are, naturally, formulated in terms of the end states u± of the shock and in terms 
of the equation of state. 

c Also, the original, one-dimensional (d = 1) derivation appears in [4]. 
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of the real part of a coefficient called 0; see (I3.38P below. That is, 

sgn Re > 

is a necessary condition for weak viscous stability. Our principal aim here is to explore, in 
a simplified setting, the effective, practical computation of 0. As pointed out by Benzoni- 
Gavage, Serre, & Zumbrun [5], the fundamental challenge for computing is the numerical 
approximation of the function y, the solution of an appropriate differential equation; see 
(I3.40p below. The determination of Re sgn has been identified as an important open prob- 
lem for physical systems due its possible role as a signal for the onset of complex behavior 
[T2] ; see also Zumbrun's more recent work on the role of the refined stability condition in 
the development of cellular instabilities for shock waves [15J. 

1.3. Plan. In §2] we describe the framework that we use as a testbed for computing y (and 
therefore 0). In particular, we restrict our attention to the case of a scalar conservation law 
with viscosity in two space dimensions. In £j3j for completeness and to make the exposition 
here mostly self contained, we recapitulate the derivation of and related Evans-function 
analysis from [5 J ; notably, because we have restricted ourselves to the scalar case, our deriva- 
tion is substantially streamlined from the general calculation presented there. In §Hwe make 
the principal contribution of this paper. We propose two distinct approaches toward com- 
puting y, and we implement these approaches in several concrete example problems. In a 
special case, we are able to find an exact solution. We use this exact solution to validate 
our numerical approximations. Finally, in §5] we actually compute and discuss the steps 
that will be required to compute using these or similar techniques in the interesting (and 
physically relevant!) system case. These calculations represent, to the best of our knowledge, 
the first ever calculations of the refined stability condition in any setting. 

2. Preliminaries 

2.1. Model. We consider the simplest possible scenario of interest, a single conservation 
law with viscosity in two space dimensions: 

d t u + d x j\u) + d x J 2 (u) = u(d 2 Xl u + d 2 X2 u) . (2.1) 

Here, u is real valued, x = (xi,X2) £ IR 2 , and t 6 1 represents time. We shall write 
a J (w) := -p-(w). The parameter v > is the viscosity; for simplicity we shall take v — \. 
Equation (12.11) is the viscous correction of the hyperbolic conservation law 

d t u + d xl f\u) + d x J 2 (u)=0. (2.2) 

Evidently, (12.21) is hyperbolic since the lxl matrix £/(u; u) = a}(u)ui + a 2 (u)u2 is always 
real for (u, uj) eKx M 2 . As is well known, such systems support discontinuous solutions 
(shocks), and we consider here the simplest possible shock solution, a planar shock wave 
connecting constant states. That is, we suppose as above that 

v/ .x /«+> x x >st 

u(x,t) = < (2.3) 

I U_ , X\ < st 

is a solution of (I2.2p . We recall that, in order for u to be a weak solution of (12.21) . u ± and s 
must satisfy the Rankine-Hugoniot condition 

s[u} = [f\u)]. (2.4) 
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2.2. Lopatinskii determinant. Applying the classical stability analysis to u amounts to 
the algebraic construction of the Lopatinskii determinant A. In general, i.e., in the case of 
( II. 2p . A is constructed from a jump term 

d 

together with bases for the stable (— ) and unstable (+) subspaces of the matrices 

A ± (\,i) = (^I + iJ2Cjdf j (u ± )j (d/Vi))" 1 . (2.5) 

For a Lax p-shock, the dimensions work out precisely since in that case 

dim E s (Aj = p - 1 , dim E U (A+) = n -p. 

However, in our setting (n = 1, p = 1), the Lopatinskii determinant consists only of the 
jump term, and we can write it down explicitly as 

A(A,0 = AM+i£[/ 2 (u)]. (2.6) 

First, observe that provided [u] ^ 0, there are no zeros of A with Re A > 0. Evidently, given 
a pair of states u± and a flux function f 2 , it is always possible to find a purely imaginary 
value of A for which A vanishes. Namely, one can simply take 



A = 



it 



The fact that there is a whole line of neutral zeros parametrized by £ is a manifestation of 
the homogeneity of A. 

Remark 1. The above calculation shows that planar shocks for (12. 2p are always weakly stable; 
see the discussion in [3]. Thus, in the current setting, the goal of computing of (3 could be 
regarded as artificial because it does not — in this case — serve to locate the transition point in 
parameter space separating the (strongly) stable and unstable regions. On the other hand, 
the simple form of A in (12.61) and the resulting abundance of weakly stable shocks makes 
this setting ideal for testing various approaches to the computation of /3. The extension of 
these ideas to a physically relevant case with n > 2 is part of our ongoing work pQ ; in §5.21 we 
indicate some features of our computations which might be useful in the setting of systems. 

2.3. Traveling-wave solutions. We now turn to the equation with viscosity. Recall, we 
have set v = 1, and we seek solutions of (12.11) of the form 

u(x, t) = u(x — st) , lim u(z) = u± . (2.7) 

2— >±00 

We write z = x — st, and we note that the traveling-wave ansatz (12 ,7p reduces the partial 
differential equation (12. ip to 

- S dz- + d-z nu)= d^- (2 ' 8) 

Integrating (12. 8p once, we find 

-s{u-u_) + f\u)-f\u_)=u', ' = A (2.9) 
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dz 



A necessary condition for the existence of a traveling-wave profile is that both states u ± be 
equilibria of the equation (12. 9p . Evidently, w_ is an equilibrium. From the requirement that 
u + be an equilibrium we recover the Rankine-Hugoniot condition (12. 4ft : 

[f\u)] = s[u]. 

Example 2 (Burgers flux). In the case that f l {u) = u 2 /2, the Rankine-Hugoniot condition 
(12. 4p reduces to 

!/ 

S = -{U+ , 

so the shock speed s is simply the average of the values of the end states. We take 

u_ = 1 , u + = — 1 
so that s = 0. Equation (12. 9p reduces in this case to 

«' = 1(^-1), (2.10) 

and a straightforward and well-known calculation shows that u(z) = — tanh(z/2). We shall 
use this exact solution of the profile equation (I2.10p to validate our numerical calculations 
below. See also Example [T2l below. 

We note that the two Lax shock conditions 

a\u + ) - s <0, a\u_)-s>0 (2.11) 

guarantee that the equilibrium at u + is a stable node and that the equilibrium at w_ is 
an unstable node. Thus, since the phase space for the nonlinear differential equation is 
one dimensional, provided that there are no equilibria between u + and w_, existence of a 
monotone profile u is immediate. Moreover, the profile will approach its limiting values u ± 
exponentially fast. 

Remark 3. By a modification of the flux f l if necessary, we may assume that s = 0. We 
make this assumption throughout the remainder of the paper. Thus, from this point forward 
the profile u is a standing wave, a function of X\ alone. 

2.4. Linearization. The first step in the stability analysis is to linearize (12. ip about the 
standing-wave solution u. We obtain 

d t v + d Xl (a\u)v) + d X2 {a\u)v) = d 2 xi v + d 2 X2 v . (2.12) 

Here, v = v(x,t) denotes the perturbation, and equation (I2.12p describes the approximate 
(linear) evolution of the perturbation v. In particular, if the linearized equation supports 
solutions v which grow in time, the solution u will be unstable. The linear equation (I2.12p 
has variable coefficients (since u is nonconstant), but these coefficients are functions of x\ 
alone. Thus, we may take the Laplace transform in t (dual variable A G C) and the Fourier 
transform in X2 (dual variable £ £ R). The transformed equation, with w denoting the 
transformed perturbation, takes the form 

\w + (a 1 ^)^)' + i£a 2 (u)w = w" - £ 2 w . (2.13) 

Here, ' denotes differentiation with respect to X\. Indeed, from this point forward, we omit 
the superfluous subscript 1 on X\\ it is the only surviving spatial variable. Thus, x will 
denote the spatial coordinate normal to the unperturbed shock front. We think of (I2.13P as 
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a family of eigenvalue problems for the collection of linear operators ££ = parametrized 
by £ G K and defined by 

:= (u/ - a 1 ^)™)' - i£a 2 (w)u> - £ 2 w . (2.14) 

We sometimes write ££{£)w = Xw as a convenient shorthand for f 1 2 . X 3 j) . Solutions of 
Jzf(£)u> = Aw which decay at ±oo with Re A > correspond to perturbations which grow in 
time since v can be recovered from w via 

v(x,t) = e xt e^w{x). 

Evidently, since u{x) — > u± as x — > ±oo, there are a pair of related, limiting operators 

& ± (£)w := w" - a l (u ± )w' - i^a 2 (u ± )w - £ 2 w . (2.15) 

Notably, for every £ G K the operators -2±(0 are constant-coefficient operators. Thus, they 
may be analyzed quite completely. This feature is an essential ingredient of the analysis. We 
also note for future reference that 

^(0)w = w" -{a 1 (u)w)' , 

and, evidently, from (I2.8P with s = 0, we find 

JSf(0)u' = 0. (2.16) 

Thus, v! is a decaying solution of (12.131) corresponding to (A, £) = (0,0). 

3. Evans function 

3.1. Evans function and low-frequency limit. The Evans function is easily constructed 
in this case. We content ourselves with a mere outline of the procedure here. For more 
details in a setting which includes ours as a special case, see [5l[T6] . First, we rewrite the 
eigenvalue problem (j£f(£) — X)w = as a first-order system of differential equations: 

W = A(x; \,£)W. (3.1) 

As above, we use ' to denote differentiation with respect to x, and we have written W = 
(w,w'y. The coefficient matrix A(x; A, £) is given by 

A ( X ' A '£)= yx + ^a 2 (u)+e + a l (u)' a l \u)) ■ (3 ' 2) 

Corresponding to the constant coefficient operators Jz? ± (£), there are constant coefficient 
first-order systems W' = A ± (A,£)W with 

A ±( A '^) = [x + ^a 2 °(u ± )+e ^K)) • (3 ' 3) 

The eigenvalues //(A, £) of A ± (A,£) are roots of 

/x 2 - fxa l (u ± ) - i£a 2 (u ± ) - £ 2 - A = . (3.4) 

Observe that as long as Re A > there can be no imaginary root of (13.41) . To see this, simply 
observe that if \x = if] for rj G R, then the left-hand side of (13. 4p has as its real part the 
expression 

which is clearly negative provided that Re A > 0. Thus, the matrices A ± have no center 
subspace on H = {Re A > 0} x ]R, and one can therefore determine the dimensions of the 
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stable and unstable subspaces of A ± by taking £ = 0. Consideration of (13. 4p then shows 
that the stable and unstable subspaces of the two matrices each have dimension one on 
H. Of particular interest are the stable subspace _E S (A + (A,£)) and the unstable subspace 
E""(A_(A, £)). The next lemma shows that there are solutions of the variable-coefficient 
problem fl3.ll) which asymptotically tend to zero (forward or backward in x) in the directions 
of these invariant subspaces. 

Lemma 4. Fix a base point P = (A ,£ ) £ {Re A > 0} x R. There exist solutions W± = 
W ± (x;\,£) of (JOD such that 

W± — >• as x -> ±oo. 

Local to P , these solutions are holomorphic in A and real analytic in £. 

Proof. The proof is a consequence of the preceding constant-coefficient analysis and the 
conjugation lemma of [10J. In particular, the conjugation lemma asserts the existence of 
a well-behaved invertible change of coordinates which maps solutions W of the constant- 
coefficient problem to solutions W of the variable-coefficient one. See [16] for details. □ 

Definition 5 (Evans function). The (local) Evans function D is defined by 

D(X, = det(W + (0; A, £), W-(0; A, 0) • (3.5) 

Our principal interest is in the polar Evans function defined for radial coordinate p G (0, oo) 

by 

D(p;X,0--=D(p\,pO. 

For ease of notation, we sometimes suppress the dependence on A and £ and simply denote 
the Evans function by D(p). 

Remark 6. The key point is to understand what happens when p — > + and when Re A = 0. 
In general, a delicate issue that arises on the boundary of H are the glancing and variable 
multiplicity sets £f and 'f . However, when n — 1, a direct computation shows that 

= 0, r = 0. 

Thus, we need not concern ourselves with these sets here. In the general case, the presence 
of glancing points, for example, prohibits the smooth extension of the stable and unstable 
subspaces of A ± to Re A = 0, and these points — a measure zero set — must be excised from 
the boundary of H in advance of the low-frequency analysis of D that follows. See [5] or [TJ] 
for more details. 

The next lemma paves the way for the fundamental low-frequency analysis of D that is 
the cornerstone of the derivation of 0. We omit the proof which is technical; details can be 
found in |16j . 

Lemma 7 (Low-frequency extension). Fix a base point P Q = (A D , £ ) £ {Re A > 0} x R. The 

Evans function D and its factors W ± = (w ± ,w' ± ) have a unique jointly analytic extension 
onto a neighborhood of (p; A, £) = (0; P a ). Moreover, the factors may be chosen so that 

w'l - {a}{u)w ± )' = 0, 

with w + (+oo) = and w_(—oo) = 0. 
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With Lemma [7] in hand, we are ready to outline the proof of the fundamental result of 
Zumbrun & Serre [16]. Importantly, it links in a rigorous way the Lopatinskh determinant 
A encoding the stability of the inviscid, ideal shock (12. 3p and the Evans function associated 
with the corresponding viscous profile ( 12. 7J . Roughly speaking, the result says that low- 
frequency/long-wave perturbations cannot distinguish between the inviscid, ideal shock u 
and the viscous profile u. 

Proposition 8 (Zumbrun & Serre, [16]). Fix (A,f) G H = {Re A > 0} x E. Then, 

D(A,0 = rA(A,0 + O(|A,e| 2 )- (3.6) 

Here, T = u'(0) is a constant measuring transversality of the connection in the traveling- 
wave problem^. The quadratic error term is uniform for (A,£) in bounded subsets of H. 
Equivalently, D(0; A,£) = and 

<y)(0;A,£) = rA(A,£). (3.7) 

Proof. We outline the principal calculation. First, recalling the discussion surrounding equa- 
tion (12.161) . we make the standard normalization at p = that to ± (0; A,£) = u' . From this, 
the assertion that -D(O) = is immediate. Second, observe that 

d p D(0) = det(d p W +} W_) + det(W+, d p W.) . (3.8) 

But then, using again the normalization at p = 0, we see that the expression for the derivative 
of the Evans function in (13.81) may be rewritten as 

d p D(0) = det(U', Y_ - Y + ) , (3.9) 

where Y ± = {y±,y' ± ) t ■— d p W± and U' = {u',u") t . In polar coordinates the eigenvalue 
equation takes the form 

pXw + (a 1 ( , u)w) / + ip^a 2 {u)w = w" — p 2 ^ 2 w . (3.10) 

Thus, differentiating with respect to p and setting p = 0, we find from (13. 1Q[) that 

Am' + (a l (u)y ± )' + i£a 2 (u)u' = y" ± . (3.11) 

We rewrite (13. lip as 

(y' ± -{a l {u)y ± )' = {Xu + ^f 2 (u))' , (3.12) 

to express both sides of the equation as perfect derivatives. We integrate the y + equation in 
( 13~T2|) from +oo to x (Recall, y+(+oo) = 0,y' + (+oo) = 0): 

y' + - a\u)y + = \{u - u + ) + i^fiu) - f 2 (u + )) . (3.13) 

Similarly, we integrate the y_ equation from — oo to x: 

y'_ - a\u)y_ = \{u - u_) + ii{f{u) - f 2 (u_)) . (3.14) 

Combining the results of (13 . 131) and (13.141) . we find that the components of Y — Y_ — Y + 
satisfy the equation 

y'-a\u)y = X[u] + ^[f 2 (u)]. (3.15) 
But, as we noted above in (12 . 16j) . u' satisfies Jif(0)u' = 0, or 

u" - a\u)u = 0. (3.16) 

^Transversality is automatic in the n = 1 case we consider in this paper. 
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Equations (I3.15P and (I3.16P show that the row operation of adding a l (u) times the first row 
to the second row simplifies the Evans determinant as follows 

d p D(0)=det(* £)=det(o ^)=TA. (3.17) 
This completes the outline of the proof. □ 

3.2. The building blocks of 0. In this section we outline the derivation of the final ingre- 
dient of 0; these calculations provide the framework for our computational approach. They 
depend on the calculations in §3.11 particularly those in the proof of Proposition [SJ The 
basic assumption is that (A ,£ ) is a zero of A with 7 Q = ReA Q = 0. In the current setting, 
this means that A Q and £ D are related via 

r = 1m\ = -£ ¥^. ( 3 - 18 ) 

Proposition 9 (Benzoni-Gavage, Serre, & Zumbrun, [5]). If(\ ,£ ) is a neutral zero of A, 
then 

/oo 
2(ir + i£ a 2 (Hv)))y(v) + ^'(v) dr) . (3.19) 
-oo 

Proof. First, we note that the assumption A = implies that (see (13.151) ) 

V ~ a l (u)y = 0. 
But, this implies that, for some constant Ci, 

y--y+ = ciu' . 

Therefore, we define 

w_ = w_ + pc\w_ , w + = w + , (3.20) 

and we note that at p = 0, we have 

w+ = w_ = u' : y_= y _ + ciu' , y+ = y+. 

We shall use this new basis for our computations. That is, we now compute with the 
alternatively defined Evans function 

D(\, = det(W + (0; A, 0, ^-(0; A, 0). (3.21) 
The advantage is that, with this new basis, 

y- = y+- (3.22) 

Then, as in Proposition [HI -D(O) = and (now also) d p D(0) = 0. We now examine the 
second derivative of D with respect to p. We observe that, by the Leibniz rule, the second 
derivative of the Evans determinant may be expanded as 

d 2 p D = det(d 2 p W + , W.) + 2 det(d p W + , d p W_) + det(W+, d 2 p W^) . (3.23) 

But, by (I3.22p . it follows immediately that 

det(8 p W +J d p W_) = det (*+ Ij J = 



y+ y. 



Thus, we may proceed in a fashion similarly as in Proposition |HJ The two remaining deter- 
minants in (I3.23P may be combined so that 

d 2 p D(0) = det(U', Z_ - Z + ) (3.24) 

where 

f) 2 W 

z ± = ■= -jj%r ■ ( 3 -25) 

We differentiate ( I3.10P twice with respect to p and we set p = 0; we find that z ± satisfy the 
equation 

2\ y ± + (a 1 ^)^)' + 2ii a\u)y ± = 2£ V , (3.26) 

or, rearranging terms, 

Jgf (0)2± = (z ± - a 1 ^)^)' = 2(A + i£ a 2 (u))y ± + 2gu' . (3.27) 

But, since y + — y_, we subtract and integrate in ( I3.27P to find that z = z_ — z + satisfies 

z' -a 1 (u)z = l, (3.28) 

where X is a constant. Next, by a row operation as in the proof of Proposition [HJ we find 
that 

figD(O) = det(U', Z_ - Z+) = det rj' = TX . (3.29) 

It remains to identify the constant X. We define ^ by = J??(0). That is, the action of 
^# on a function g is given by 

J?g{x) = g'(x) - a l (u(x))g(x) . 
Thus, we may rewrite (I3.28j) as 

X = Jtfz_(0) -J?z+(0) . (3.30) 
But, by the fundamental theorem of calculus, 

JS?(0)z_(x) da; = JlzJp) - ^z_(-oo) , (3.31) 



and 

/" + OO 

Sf(0)z+(x) dx = JZz + {+oo) - ^z + (0) . (3.32) 



Thus, using ( I3.3ip and (13.32)1 . we are finally in a position to identify the constant X. We see, 
starting with the expression in ( I3.30p . that 

X= Jtz.{$S) -J?z+(0) (3.33) 

p+oo 

5?(0)z_(x)dx + ^f(0)z + (x)dx + ^z_(-oo) - Jtz+(+oo) (3.34) 



oo 

+ 00 



Sf(0)z ± (x) dx + Ji~z_{-oc) - ^#z+(+oo) (3.35) 



oo 
+oo 



2(A D + \i a\u{x)))y{x) + 2^u'(x) dx . (3.36) 
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We note that the equality in ( I3.35P follows from the fact that j£f(0)z + = J2?(0)iL (see 
(I3.27P ). and the boundary terms vanish because z± tends to zero exponentially as x — > ±00. 
In conclusion, we have shown that 

/+00 
2(A D + i£ a 2 {u{x)))y{x) + 2gu'(x) dx . (3.37) 
-00 

□ 

Propositions[H]and[H]provide the fundamental ingredients for the derivation and calculation 
of /3. The derivation is based on the relationship between the zero sets of A and D expressed 
in Proposition [HJ The main analytic tool is the implicit function theorem. We omit the 
details, which may be found in [16], and simply note that f3 is given by 

d 2 D\ r d 2 D\~ l 



p= W)\m) ■ (3 ' 38) 

But, Proposition [8] implies that 

3 2 D r ^A 



dpdX d\ 

and — given the explicit form of A in ( 12. 6ft — we see immediately that 

dA 

Thus, from Proposition [9], we find that 

= (3.39) 

From ( I3.39P we see immediately that the computation of (5 requires one to compute X. 
Evidently, from ( I3.37p . to compute I one needs to know both the profile u and the function 
y solving the differential equation (see (I3.12p ) 

(y' - a l (u)y)' = (A D + i£ a 2 (u(x)))u' . (3.40) 

In the next section we propose two methods for the practical approximation of y. 

Remark 10. Briefly, the sign of Re /3 detects whether or not the zero level set of D curls 
into the unstable half plane. Heuristically, it has a physical interpretation as an "effective 
viscosity" coefficient for transversely propagating deformations of the front. A detailed 
discussion of this point can be found in |16j . 



4. Computing y 

As pointed out in [5], the principal task in finding /3 is to compute y; the other elements 
required to compute /3 are the building blocks of A — in general these are the eigenvalues 
and eigenvectors! of A±(\,£) — and the profile u. Benzoni-Gavage, Serre, & Zumbrun [5] 
have proposed a two-step method for finding /3. That is, the equation for y is a linear 
equation whose coefficients depend on the profile u. Thus, they propose to first solve the 
profile equation ( 12. 9p . For gas dynamics, techniques for doing this are well known; see, 
e.g., [SI. Then, they describe how to transform the equation for y into one that fits into a 



e In the general case, one also uses these eigenvalues and eigenvectors to compute a boundary term B; this 
term does not appear in our setting. 
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standard numerical framework for approximating the solution of a two-point boundary value 
problem. Their method has never, to our knowledge, been implemented. Below, we propose 
two alternative approaches that, at least in the present context (n = 1), work well. We are 
currently investigating all three approaches in the case of gas dynamics for which n > 1. 

4.1. Integrating factor. This method exploits the linear structure of (I3.40p . The equation 
for y is 

(y' - a\u)y)' = (ir + i£ a 2 {u(x)))u' . (4.1) 
We integrate both sides of (14. ip . and we obtain 

y' - a\u)y = ir (u - u.) + i£o(/ 2 (fi) - / 2 «>) . (4.2) 
We write y in terms of its real and imaginary parts as 

y — w + iv . 

Then, the relevant system of equations becomes 

w' = a 1 (u)w , (4.3a) 

v > = a \u)v + r (u - u_) + Uf(u) - f(u_)) . (4.3b) 

We suppose that u is known, and we write 

F(X) = T {U(X) -U_)+ Uf{u{x)) - f{u_)) . 

Now, suppose that M is defined by 



M(x) := exp f- J 



a 1 («(z))dz . (4.4) 



Then, a simple calculation shows that M'{x) = — a l {u{x))M{x) 1 whence the equation for 
w can be rewritten as a perfect derivative (Mw)' = 0. Integrating from to x and using 
w(0) = A, we find immediately that 



w(x) = Aexpyj a (u(z))dzj . (4.5) 

Thus, w is given explicitly in terms of the profile u. Similarly, we apply the integrating factor 
M to the equation for v. We write v (0) = B, and we see that 

v (x) = e tf al W*» d * | B + J\-^ al ^ d ^{z) &z 

e /; a> (u(z)) dz | B + I* e - S» ^ (U(V)) d V ^ ^ _ u _ J + ^ (/ 2 ( - ^ j _ f 2 ^_ j j ^ 

# L6) 

Equations (14. 5 p and (14.61) show that we may write y = w + iv explicitly in terms of the profile 
u. Thus, rather than solve a differential equation for y, in this case we may simply approx- 
imate u (as described above, a well understood problem) and then use that approximation 
to compute the integrals in (I4.5p and (14. 6p . Indeed, as Example [T21 below shows, in some 
cases it is possible to compute these integrals exactly. 
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Remark 11. From [5], we note that in the computation of 0, the two apparently free param- 
eters A and B above should be chosen so that y satisfies an orthogonality condition at the 
origin. In the current setting, this reduces to the requirement that 

A = 0, B = 0. 

solution). In the case of the Burgers flux = u 2 /2, we have seen 

= 1, then s = 0. We have also seen that in this case the profile u is 

u(x) = — tanh(x/2) . 

Suppose that f 2 {u) = u 2 . We take r D = and £ = 1 so that A(ir ,£ ) = 0. Then, system 
of equations (14.31) reduces to 

w' = uw , v' = uv + u 2 — 1 . (4.7) 

That is, we need to solve 

w' = — tanh(x/2)w , v' = — tanh(x/2)t> + tanh 2 (a;/2) — 1 . 
By direct computation, with A = and B = 0, we find immediately from (14. 5p . (14. 6 p that 

w(x) = , v (x) = — x sech 2 (a;/2) . (4.8) 
This solution is plotted in Figure El 

Remark 13. The formulae in [5] are derived in the case of strictly parabolic viscosity. That 
is, the matrices B^ k in (ll.3p are assumed to satisfy 

Rea(j2 ^kB jk (u(-)) \ > for all £ e R d \ {0} . (4.9) 
\j,k=i J 

However, the important physical case of gas dynamics features only a partially parabolic or 
"real" viscosity. For example, consider the equations of isentropic gas dynamics 

d t p + V- (pu) = 0, (4.10a) 
d t (pu) + V • (pu eg) u) + Vp = pAu + (p + rj) graddivu , (4.10b) 

where p and rj are the first ("dynamic") and second viscosity coefficients (assumed here to 
be constant). The lack of second-order terms on the right-hand side of the conservation 
of mass equation (I4.10al) prevents the system (I4.10p from satisfying (14.91) ; nonetheless, it 
is clear that the results of [5] extend in a natural way to systems with "real" or partially 
parabolic viscosity, such as the equations of gas dynamics. It is a rather tedious exercise to 
derive the equation for y in that setting. Remarkably, our preliminary calculations for the 
system (14. lUp show that the equation for y can be written as a linear diagonal system [TJ; 
this suggests that the above approach based on integrating factors might be applicable to 
the corresponding calculation for (14.101) . 



Example 12 (Exact 
that if u + = —1, u_ 
given by 



4.2. Coupled formulation. 



13 



4.2.1. Description. Rather than solve the problem in two steps (first u and then y), we 
consider now the problem of solving the coupled system for u and y = w + iv. Thus, we 
consider the autonomous system 

u' = f\u) - f l (uj , (4.11a) 

w' = a\u)w, (4.11b) 

v > = a \u)v + r (u - u_) + Uf(u) - f 2 {u_)) . (4.11c) 

Our motivation for pursuing this approach is that, if the numerical methods involved do not 
take explicit advantage of the linear structure of the y equation, for modestly sized systems 
it involves no extra work to solve the coupled system in a single step. On the other hand, 
the challenge in this case is to construct a suitable guess for the solver. However, once a 
suitable guess is found, the calculation can take advantage of continuation. For example, 
if u_ is being moved along the Hugoniot curve, one can use the previously found solution as 
the initial guess for the next value of u_; see Example [TBI Based on our experiments, this 
method works quite well. 

We write (14.111) as U' = F(U), where U = {u,w, z). The desired solution is a heteroclinic 
orbit in the phase space R 3 which connects the equilibria U± = (u ± ,0,0). The linearization 
is straightforward; we see immediately that 

/ a}(u) \ 

dF(U) = a\u) . (4.12) 

V§^ + r + £oa 2 (n) a\u)) 

Evidently, the eigenvalues of dF(U±) are, with multiplicity three, simply a 1 (u±). From ( 12. lip 
this gives a connection from a stable node to an unstable node in M 3 . We expect then to 
introduce two parameters to completely parametrize the solutions; as noted in Remark [TTj 
an orthogonality condition will select a particular solution for the computation of /3. 

4.2.2. Numerical Implementation. We truncate the problem to the computational domain 
[—L,L\. Thus, the system becomes 

U' = F(U), xe[-L,L]. (4.13) 

At this point we expect to incorporate projective boundary conditions at ±L. Next, for 
convenience, we double the variables and rescale the problem to the unit interval. Thus we 
consider the problem 

U' r = LF(U r ) , U' t = -LF(U e ) , x G [0, 1] , 

where U r (x) = U(Lx) and Ue(x) = U(—Lx). We also implement the classical phase condition 
to remove the translational invariance from the problem. Folding over the solution makes 
it simple to include this as a boundary condition at x = 0. Thus, in this case, we use 
one phase condition, three matching conditions, and two free parameters to determine the 
solution. The penultimate step is to generate an approximate solution, or a guess. In this 
case, we generate our guess by solving the corresponding initial-value problem. In practice, 
since this is a sink-source connection, it is not too difficult to generate a good guess by this 
method. Once a guess of sufficient quality is found; solutions for nearby parameter values 
can be found easily by continuation. Finally, we solve the boundary-value problem using 
MatLab's routine bvp5c; this is a code that implements the four-stage Lobatto Ilia formula; 
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this is a collocation formula. The collocation polynomial provides a C 1 -continuous solution 
that is fifth-order accurate uniformly in the computational domain [0, 1]. Sample solutions 
using this procedure are plotted in Figure [1] and Figure [2j 

Remark 14. A natural question concerns the determination of the smallest size of L that 
guarantees that the numerical approximation fully resolves the features of the true problem 
on R. Here, for the model problem that we consider, we are content to do so in an ad hoc way. 
Since our ultimate interest is in sgn Re (3, we merely verify that the value of this quantity is 
stable as L is increased; we take this as evidence that the computation is sufficiently resolved. 
See Table [Q 

4.2.3. Examples. 

Example 15 = u 2 /2, f 2 (u) = u 2 ). This is the same as Example [T2"l and so there is 

an exact solution available. The system takes the form 

u' = ^(u 2 -l), (4.14a) 

w' = uw , (4.14b) 
v' = uv + u 2 — 1 . (4.14c) 

The approximate solution computed using the method described in Section 14.21 with L = 20 
is plotted against the exact solution in Figure [2j We note that the 2-norm error between the 
approximate solution and the exact solution on their common domain is well-controlled by 
the built-in error control features of bvp5c. 

Example 16 (Moving end state). In this example, we take as before f l {u) = u 2 /2, and we 
take 

f 2 (u) = sin(47ru) 

for the transverse flux f 2 . Our aim in this example is to mimic the kind of calculation that 
one would do in practice, searching for the point in parameter space (e.g., the value of u_) 
at which sgn Re (5 changes sign. To that end, we systematically increase the value of u_ and 
recompute the u and y for each new value of the end state. This is precisely the kind of 
computation that is well suited for continuation. The results are plotted in Figure [TJ 

4.3. Comparison of two methods. In the context of the scalar conservation law (12. ip . 
both methods work well. Our experiments show that some care should be exercised in 
approximating y via the integrating factor method; in particular, the overall quality of the 
computation depends on the approximation of the integrals in (14. 6p . Our implementation 
uses Simpson's rule to approximate these integrals. It is worth noting that the profile u 
is computed independently in this first method; for example, the calculation does not take 
any account of the transverse flux f 2 . By way of comparison, error control in the coupled 
formulation is "automatic" since we use built-in convergence tolerances in the package bvp5c 
to control the quality of the approximation across the entire computational domain; That is, 
u and y are treated on the same footing, and this method takes account of the entire structure 
of the problem at each stage of the iteration. In Figure [2] we plot the exact solution from 
Example [T2] against the approximate solutions obtained by the two methods. The caption 
of that figure records the error between the approximate solutions and the exact solution in 

15 



V 



(d) ' " ' •(e) - ' • (f) ' " ' " 

Figure 1. The effect of moving the left state. Black: profile u, Blue: w = 

Rey, Red: v = Imy: (a) u_ = 1. (b) u_ = 1.1. (c) u_ = 1.2. (d) u_ = 1.3. 
(e) u_ = 1.4. (f) u_ = 1.5. 



the 2-norm, given by 

^ .v 

\ x h 



J>J. (4.15) 

\ 3=1 



On the other hand, a distinction of note between the methods is in the complexity of 
the nonlinear two-point boundary-value problem that must be solved. Our technique for 
solving these problems by collocation hinges on finding or computing a suitable initial guess. 
This guess is used as the seed in an iterative solution of the nonlinear equations for the 
coefficients of the collocation polynomial. In the first method, based on the integrating 
factor, a relatively simple boundary-value problem needs to be solved; one expects that it 
is, generally, much easier to generate a good initial guess for this problem than to find a 
similarly good initial guess for the coupled formulation of the problem. 

5. Conclusion 

5.1. Calculating 0. Finally, with our approximations of u and y in hand, we proceed to 
compute p. Examining (I3.39p . we observe that the formula for j3 may be rewritten as 

rx 



P 



FA, 

[u]- 1 I 2(ir D + i£ a 2 (u))(w + iv) + 2gu' dx 



COD 



DC 



2 ' L 



(ir + it a 2 (u))(™ + iv) + Cifiu) - f\u_) - s(u - u_)) dx . (5.1) 

U\ J-L 

Thus, from §4.11 and §4.21 we have two distinct ways of obtaining approximations for u, w, 
and v on the computational domain [— L, L] in order to approximate the integral in (15. lft . 
In Table [T] we compare the values of /3 computed by each of the methods for the problem in 
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FIGURE 2. The computed solutions (autonomous = stars, integrating factor 
= circles) plotted against the exact solution (solid lines) from Example [121 
Here, L = 20, and ||u auto - u exact || 2 =1.0470e-07, \\w auto - u> exac t||2=0, and 
ll^auto -fexact||2= 4.42128e-07. Similarly, ||u int -u exact || 2 =1.45990e-06, \\w int - 
, u ; e X act||2=0, and 1 1 Ui nt — f ex act||2= 2.79917e-04. The error in the autonomous 
formulation is controlled directly by the convergence criteria of bvp5c while 
the error in the IF formulation depends on the accuracy of the computed profile 
u and the quality of the approximation of the integrals in (I4.6p . 



L 




10 


20 


30 


/^auto 


9.9918 


+ O.OOOOi 


10.0000 + O.OOOOi 


10.0000 + O.OOOOi 


Ant 


9.9919 


+ O.OOOOi 


10.0001 + O.OOOOi 


10.0001 + O.OOOOi 



Table 1. The computation of (3 for Example [T2l /3 ex act = 10. 

Example [T2l In each case, we approximate the integral in ( 15.1 ft by the trapezoid rule, and 
we denote the different values obtained by /3 auto and f3 int . 

5.2. Discussion. We have considered Zumbrun & Serre's refined stability condition [16] 
in the simplest possible setting, and we have proposed and implemented two methods for 
practically computing the condition. Both methods work well in the present setting, and we 
believe that the approaches for computing y from §4] are potentially useful for interesting 
physical versions of this problem. Indeed, our principal interest is to use the model problem 
here as a stepping stone towards the analogous problem for a physical system (n > 1) 
which possesses open set of weakly stable waves. For isentropic gas dynamics in two space 
dimensions ()4.10p . an example of such a family has been given by Majda [9]. In that case 
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n = 3, and the fundamental new issue that arises is the presence of slow modes, see [HE]. 
For example, in the case of a Lax 1-shock, there are two outgoing characteristic directions, 
and this alters the calculation from the very beginning. For example, the computation of 
A is more involved (although it is still computable in closed form In this case the 

Lopatinskii determinant takes the form 

A(A,0 = det(r 2 + (A,0,r 3 + (A,0, A[f/] + ^[f(U)}) . (5.2) 

where the vectors {rj (A, £), (A, £)} are a basis for the unstable subspace of A+ where 

A ± (X,0 = (xdF°(U ± ) + i£dF 2 (U ± ))(dF 1 (U ± ))- 1 , (5.3) 

with U = (p, ■u 1 ,M 2 ) t and 

/ P \ ( PUI \ ( PU2 \ 

F°(U) = \pu x , F\U) = pu\+p{p) , F\U) = pu lU2 . (5.4) 
\pu 2 J \ pu x u 2 J \pul+p(p)J 

In the fluxes , the form of the pressure p is prescribed as the equation of state — the 
constitutive relation that specifies the nature of the gas. Mathematical treatments of gas 
dynamics frequently take a "7-law" gas^ 

p{p) = a p 7 , (5.5) 

where oq is a positive constant and 7 > 1. For such a pressure law, there are no weakly 
stable shocks [11]. However, thermodynamically admissible perturbations of a 7-law pressure 
function can open up regions of weak stability [HE] • Additionally, the presence of slow modes 
introduces an additional boundary term £>, as described in [5], into the formulation of (3. 
The term B is associated with the absence of a spectral gap for the linearized operator J?f (0), 
and B is constructed from the right and left eigenvectors of A ± in (15. 3 p ; in the case of (I4.10p , 
these are known explicitly. Finally, in the physical cases, low-frequency behavior on OH can 
be complicated due to the presence of nonempty glancing and variable multiplicity sets <S 
and "V . In the case of (14.101) "V is empty and £f can be computed explicitly. The application 
of the ideas in §Hto this problem is a topic of our current investigation. 
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